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The mathematical model of a rotating electrohydrodynamic flow in a thin suspended liquid film 
is proposed and studied. The motion is driven by the given difference of potentials in one direction 
and constant external electrical field iSout in another direction in the plane of a film. To derive 
the model we employ the spatial averaging over the normal coordinate to a film that leads to the 
average Reynolds stress that is proportional to |i3out|^. This stress generates tangential velocity 
in the vicinity of the edges of a film that, in turn, causes the rotational motion of a liquid. The 
proposed model is aimed to explain the experimental observations of the liquid film motor P, 0]. 

PACS numbers: 47.32.Ef, 68.15.+C, 47.57.jd 
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I. INTRODUCTION 



The paper is devoted to the theoretical study of the 
motions in a thin suspended hquid fihii. The motions 
are driven by the constant external electric field that is 
applied at the edges of a film. We show that this field 
produces the averaged rotating motion of the hquid in the 
plane of a film. Our aim is to explain the rotating flow 
observed in an electrolyze planar water cell placed inside 
a plane capacitor [l], Q . The rotating motion of a fluid 
as a whole caused by the action of a constant electrical 
field is so unusual that the authors of P, call this effect 
a liquid film motor, emphasizing that it represents a new 
type of engine. They also proposed that it could be ex- 
plained by the changing of orientation of water molecular 
dipoles caused by a strong electric field. Simultaneously 
they denied the possibility of the generating of such a 
flow by the edge effects. In contrast, we show that the 
jump of an electric field across a water-dielectric inter- 
face produces the tangential velocity of a liquid that can 
maintain a steady rotating flow in the whole film. In 
other words, we demonstrate that one docs not need to 
use a heuristic idea about the switching of the molecular 
orientations: this phenomenon can be explained with the 
use of classical tools only. In our model the rotating mo- 
tion in a film is explained by the electro-kinetic effects at 
its edges. According to our theory the ratio between the 
spatial scales of a flow domain plays a crucial role. Only 
for thin fllms the classical edge effects can generate the 
rotation; in contrast in the flow domains with all spatial 
scales of the same order this effect will be absent. Nat- 
urally, our flnal model is two-dimensional (plane); how- 
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ever the tangential velocity at the boundary is actually 
caused by Reynolds stresses that appear after the averag- 
ing over a fllm thickness of an original three-dimensional 
flow. The resulting tangential velocity has order 0{h'^), 
where 2h is a fllm thickness. An intense electrohydro- 
dynamic (EHD) rotating flow takes place only in the re- 
stricted domain of governing parameters. In particular, 
such a flow can exist for the fllms of moderate thickness 
(for example, the tangential velocity ^ 1 cm/sec appears 
for the following parameters: the strength of the capaci- 
tor electric fleld ^ 30kV/m, the difference of electrolysis 
potentials ^ 20 V, the fllm thickness ^ 0.1-0.3 cm, and 
the fllm surface size ^ 1 cm) , but it can not exist for very 
thin films. An important general result of our paper is 
the demonstration of the fact that the classical effects 
(such as the electrokinetic phenomena), that are small 
in ordinary conditions, can play the key part in micro- 
scales. Such revaluation of the classical effects may be 
important for the developments of microfluidics and for 
the creation of microdevices. Here one have take into 
account that in this paper we both present systematical 
theoretical results and show their good agreement with 
the experiments of P, • The detailed discussion of our 
results is given in Sect. lVIll 

In the mathematical modelling we essentially use the 
results H, 0, H, 0, 0, [ll that contain the analytical and 
numerical studies of EHD flows with the gradient of con- 
ductivity, the method of depth-average, and an effective 
asymptotic procedure for the EHD equations of multi- 
cornponent mixtures. In the averaged equations derived 
in 0, H, @, 01 one can sec the terms correspondent 
to Taylor-Aris dispersion and to Reynolds stresses; how- 
ever Reynolds stresses are neglected since they are small 
for the chosen intervals of parameters and negligible for 
the studied phenomena. It is also important that in 
one can flnd the comparison between the results for 
the mathematical models of the different levels of ap- 
proximation. The papers, closely related to our studies, 
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@,[i3,[ni,[ll,[ll consider EHD flows in thin liquid films 
or in liquid layers with interfaces; paper [iJl describes 
the appearance of vortex rings due to reactions near an 
electrode; and presents a rotating EHD flow in a 
smectic medium. The role of interface boundary condi- 
tions in the electrohydrodynamics (EHD) is well-known 
from the classical papers [l^, [13, [3- The paper [l^ 
shows that an electrical double layer (EDL) can allow 
a slip in the boundary conditions between a liquid and 
a solid. The papers ^ iH, [H, [H, [13, ^ are devoted 
to the influence of the inhomogeneous electrical charge of 
microchannel boundaries on EHD flows. Other closely re- 
lated papers [H, [S^, [H, [2^, [sO] consider various theories 
of the EDL, including so-called extremal regimes. The 
survey of modern EDL-theories can be found in [3l|. 

The interest in various flows of micro- and nano-scales 
has increased greatly during the last few years. For exam- 
ple, the main parts of the surveys arc devoted to 
EHD processes in microchannels; .31, .38, 
deal with the injection of a fluid and other flows in mi- 
crochannels. This interest is strongly stimulated by the 
creation of the microfabricatcd fluid devices for the ef- 
fective separation or the micromixing of multicomponcnt 
mixtures 0, [HI, [13, [H, [13, [H , the electro-micro-pumps 
[1^ , etc. These new techniques are known as parts of the 
Lab-on-a-Chip technology. The mathematical models in 
this research help to understand and to describe 
micro-processes, to develop experimental methods, and 
to construct microchips. 



We assume that the electric fleld is potential, the grav- 
ity and the surface tension are absent; the dielectric per- 
mittivity e = const that leads to the absence of the pon- 
dermotive force (l/2)Ve(V(/3)^ = 0. The dimensionless 
system of governing equations that describes EHD flows 
of a multicomponcnt fluid (for example, water with the 
ions H+, 0H-) is: 

du 

S^— = -5'^\/oP + S^i^Aou + ud^.u - gVoV?, (1) 
at 



at 



divo u + dzW = 0, 



(3) 



e{S'^Aoip + dzz'p) = -S'^q, q = ^ekCk, (4) 



i5'^+(5'divoife + a,/fe =0, (5) 



ik = -£'fc(VoCfc + efc7CfcVo'^), (6) 



h = -DkidzCk + CkjCkdzip), 



T = {ai3,a23,0)^HdzU + S^Wow), (7) 



II. BASIC EQUATIONS 

The rectangular thin liquid film with the fixed plane 
free surfaces z = ±h is considered in Cartesian coordi- 
nates {x,y,z) (Fig.[T|). The electrical field can be conve- 
niently split into two parts. The first one is due to the 
constant electric potentials f = and ip = ip o on the 
boundaries a; = and x = X, so the constant difference 
of potentials is applied in the direction a;. The second 
part is the constant external electric field £^out that is 
prescribed at the boundaries y = and y ^ Y. The 
vector Eout hes in the plane z = const, a is the angle 
between this vector and y-axis. 




FIG. 1: A thin film. Dimensionless variables are given in 
parenthesises. 



Here v = (it, w) and u = (u, v) is velocity and its (cc, y)- 
projection, p is pressure; q is the density of molar charge; 
(p is electric potential; is the molar concentration for 
the fc-th component of mixture; ik and Ik are the planar 
and transversal density fluxes for the concentrations; ly 
is kinematic viscosity; Dk is the diffusivity for the com- 
ponents of mixture; Ck are the electric charges of compo- 
nents (in the units of electron's charge) ; e is the solution 
permittivity; the parameter 7 characterises the ratio be- 
tween the transports of concentrations by the electrical 
field and by diffusion; 2S is a dimensionless film thick- 
ness; r is the tangential stress vector that is expressed 
via the components ctis, a23 of a viscous stress tensor. 

On the film boundaries z = ±1 we accept: the no-leak 
condition for velocity 

the stress-free condition, that with the use of ([S]) is 

rL=±i=f^(5.M + <^'Vow;)|^^±,= t^a,«|^^^^==0, (9) 
the no-leak conditions for concentrations 

4L^±i=0 (10) 
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and the vanishing of the normal electrical current 



±1 



0. 



(11) 



We use the governing equations ([T|)-([7l) for the deriving of 
the averaged model in Sect. lIIII and in AppendixIXl In our 
averaging procedure we use only the boundary conditions 
([51)- pT|) . The other boundary conditions (defined for the 
averaged equations) are given in Sect. lIIII 

For the introducing of dimensionless variables we use 
the following characteristic values of parameters: 



[z] = h, [t] = T, [u, v] 



a 



[cfc]=C, [E]=8, [p]^£a, M=FC, 7 



F£a 
'RT' 



[p] = FCSaS'^, 



pa 



FC£5^ ' 



-T- (12) 



Here a is the characteristic length in the plane of the 
film; h and 5 are the dimensional and dimensionless half- 
thickness of the film; p is the density of a liquid; T, P, 
C are characteristic time, pressure, and molar concentra- 
tion; FC is characteristic charge density; F is Faraday 
constant; R is the universal gas constant; T is the ab- 
solute temperature of solution; a£ is the characteristic 
difference of electric potentials in the x-direction. The 
dimensional values of kinematic viscosity v* , diffusion co- 
efficients D^, and the dielectric permittivity e* are linked 
to their dimensionless counterparts as: 



Dk 



Dir 



e = 



e*£ 
oFC' 



(13) 



The use of dimensionless parameters (viscosity, diffu- 
sivity, etc.) instead of conventional scaling numbers 
(Reynolds number, Peclet number, etc.) is more con- 
venient for our purposes since they allow us to see what 
physical effects participate into a certain process. The 
connections between the introduced parameters and the 
scaling numbers are apparent: 



Rc 



1 



Fez, 



1 



III. THE AVERAGING ACROSS A FILM 

The main part of the employed averaging procedure is 
the same as in 0, [4, 5, 6, 7J. The operation of averaging 
is defined as: 



fix, y,t) = ^ I fix, y, z, t) dz, f = f- f. (14) 



We decompose the solution of ([T|)- pT|) into the series 
{u, w,p, q, cu, ^} =E{"" J'" 9" ^"i-^'™ = 



Tn=0 



{u ,11) ,p ,q ,Ck,tp \6 . 



(15) 



m=0 



The averaging of the governing equations (HJ-®, which 
takes into account the boundary conditions ([8|)- (fTT|l and 
the decomposition in small parameter 5^ , first yields q = 
t + Oi5^), ^ = ^0 + 0(<52), cfe = c° + 0(<52), ^0 =^0, 
= 0, = and then leads to the expressions for 

The averaged plane equations which keep the 



w 



~i 



terms 0((5^) are (for the details see Appendix E| 

dpi 
dt 



-5^^ + /3Vo(«7 ®U) = -S^VoP + S^i^Aqu - i^J7,(16) 



divo M = 0, 



eAoip^-q, q = ^ekCk, 



(17) 



(18) 



dpCk 
dt 



akS^ divo(C/(J7 • Vocfc)) + divo ik = 0, (19) 



ik = -£'fc(VoCfc -I- CkjCk'^ov) 



do 
dt 



dt + u- \/o, vU = qS^of, 



(20) 



^ 45' 



ak = 



94:5Dk ' 



where (C/ ^ U) denotes a tensorial product. We em- 
phasize that after this averaging S must be treated as a 
regular independent parameter of the problem, together 
with ly, e, Dk, etc. 

For the equations P^ - (|^D|) we prescribe the boundary 
conditions for the averaged fields u, Ck, ^ on the side 
boundaries x = Q,X and y = Q,Y (Fig.[T]). 

The boundaries y = and y = Y represent the inter- 
faces between two dielectric materials: the liquid with the 
dielectric permittivity e and the outside medium with the 
dielectric permittivity Eout; these boundaries are insula- 
tors (not electrodes), hence we prescribe the continuity 
of the normal components for electrical induction [13, l48j 

e(n • Vo<^) = eout(ri • £?out), y = Q,Y, 

where n is the unit normal vector to the boundary. Since 
the vector £^out lies in the plane z = const and has the 
angle a with j/-axis we have 



dip 
dn 



I^Joutlcosa, (21) 
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where the sign '— ' corresponds to the boundary y = 
(Fig.[T]) [5^. The conditions of zero concentration fluxes 
at y = 0, y are 



IV. THE FLOWS NEAR BOUNDARIES 



ik-n = 0, y = 0,Y. 



(22) 



The fixed difference between the electric potentials at 
X = 0, X = X is given as 



ip = 0, X = 0; ip ~ ipQ, X ~ X. 



(23) 



For all edge boundaries x = 0, X and y ~ 0,Y we require 
the no-leak of a liquid 



la:— 0, A ' ly— 0,1 



and the conditions 



u ■ t\ ^ ^— —TZWop ■ t\ „ ^, 

I y— 0, r ^~ ly— 0,r ' 



=0, A' ' 



(24) 



(25) 



(26) 



where t is a unit tangent vector to the boundary, TZ is 
the coefficient defined in Sect. lIVI By virtue of ((23|) the 
boundary conditions for x ~ 0,X take the form of 
no-slip condition 



(27) 



The prescription of the tangential velocity ([25|) at the 
boundaries y = 0,Y is justified in Sect.[lV] where we 
derive it and show that TZ ^ Eq. It is derived from 
the equation (flB)) that contains the averaged Reynolds 
stresses 



pVniU ®U)= p{U -VqU + UAiYaU), (28) 

which define TZ for certain intervals of parameters v, S, 
e, etc. 

The derivation of p^ -([20 |) is given in Appendix [Aj 
here we mention only that the boundary conditions ([5])- 
pT|) at z = ±1 play a central part in this derivation. 
It is also well-known that if we use only spatial averag- 
ing it does not allow us to produce the closed systems of 
equations; for its closure one has to employ some addi- 
tional hypothesis. As such a hypothesis we propose the 
condition = that is physically natural and accepted 
without any mathematical justification. The equations 
similar to (HHl-dlDD have been obtained in 0, |3, M 0, 
(and other papers cited there) devoted to the studies of 
EHD flows with the spatially nonuniform conductivity in 
the microchannels with solid boundaries. These papers 
also contain the decomposition into the power series and 
even the term similar to ([28)) . However the key difference 
with ourpajper is: for the physical parameters considered 
in 0, 0, H, @, 01 this term is small, so it is naturally ne- 
glected. 



The problem p6l) - ([28)) can be split into the sequence 
of two problems: (i) the calculation of TZ in ([25|l . and (ii) 
the finding of the averaged velocity field u and the poten- 
tial Tp. In order to evaluate TZ we assume that the mixture 
is clcctroneutral everywhere except the vicinities of the 
boundaries y = 0, y. In these vicinities we build the 
boundary-layer solution that leads to a good estimation 
for TZ. The detailed studies of the related double layers 
(the GouyChapman layer, the Stern layer e tc.) can be 
found in 0, m, US, S, S, 113, H [13, [111, where 
nonlinear and steric effects are taken into account along 
with the linear electrokinetic effects. From the math- 
ematical viewpoint different EDL theories are aimed to 
formulate and to justify different boundary conditions for 
the related boundary layers. The main question is how 
to choose the mutual positions of a physical boundary 
and an interface between the regions with positive and 
negative charges. 

Let us consider the vicinity of the boundary y = 
(the case y — Y is similar) and look for a steady solution 
of the problem (fT6|) -([24 |) neglecting in (fT9|) the terms 
afe52[7(j7 . Vgcfe) 

ll=(«(y),0), Ck = ck{y). V = Hv) + Ex, (29) 

where E is the constant tangential component of the elec- 
trical field in the vicinity of y = 0. The integration of 
with the boundary conditions yields 



Ck{y)=CBke q{y) =^ekCBke 



(30) 



where cbu are the constants representing concentrations 
for the equilibrium Boltzmann distributions. We restrict 
ourselves with the case when the mixture is electroneu- 
tral; only two kinds of ions arc present (for example 
and OH" for water), and the equilibrium Boltzmann dis- 
tribution is valid: 



CBl 



dct 

CB2 = CS, 



Zl 



1, 



Z2 



-1. 



(31) 



The Poisson-Boltzmann equation p8|) takes the form 



X^dyyO^sinYiB, e{y)^-f<^{y). 
where A_d (or A) are (or relative) Debye's length: 

e*RT 



2jCB 



< 1, 



>2 - 



2c*„F2 ■ 



(32) 



(33) 



In the vicinity of y = the boundary layer variable is 
introduced as y — Xr] (similarly, at y = y the change of 
variable is y = Y + Xrj). In the more precise terms the 
considered boundary-layer solution represents so called 
'penetrating boundary layer' [i^. In this case the origi- 
nal equations and the boundary-layer equations coincide, 
and A <C 1 is not required for the obtaining of the solu- 
tion. Instead of the looking for the boundary-layer solu- 
tion decaying at infinity one can use the symmetry with 
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respect to the middle of the domain (y = F/2); the result 
will remain the same. Nevertheless our further consider- 
ation follows the path that is more transparent from the 
physical viewpoint. The equation (|32p takes form 



dtn^iv) = sinhe'(77). 



(34) 



Its integration with the boundary condition (|2T|) for $(0) 
yields 



9(0) = -00, 7$(0) = -^ 



(35) 



where 



9o = In f 1 + + Ex/2 + eA , = —E^, 



The expression for ^(Y) is similar to ([25 

6(0) = 00, 7$(y) = 00, 



where different signs correspond to j/ = and y = Y. 
With a sufficient precision the value of TZ^ at E ^ 1.5 
can be replaced with 



± 



E 1.5. 



(39) 



In order to avoid misunderstanding we should mention 
that the calculated value of TZ (^5]) represents only a 
rough estimation; to obtain it we accept that the equa- 
tion is steady and neglect the Taylor-Aris disper- 
sion. Moreover, in ([29|) we assume that E = const on 
the boundaries y = 0,Y that is not true, later on we 
consider E = E{x) (sec p5|) '). In spite of these simpli- 
fying assumptions, the results of this section show that 
Reynolds stresses p\/o{U(E)U) for certain parameters can 
crucially participate to the generation of the tangential 
velocity (of order 0{Eq)) at the side boundary of a film. 



The boundaries y = and y = Y represent the different 
plates of the capacitor, therefore the opposite signs of the 
potential $ are apparent. The calculation of /3Vo(C/(3)C/) 
with the use of (I2H1), (EOl), (EH) yields 



/3Vo(?7 ®U) = ^ (Edy {tdy<^) , dy [qdy^)' 



(36) 



where the righthand side (written in components) allows 
us to integrate the equation p6|) with the additional con- 
dition u(oo) = 0. This condition means that the fiow 
arising near the boundary must decay at large distances, 
i.e. that the distributions of the horizontal component u 
and potential $ are of a boundary-layer type (for more 
details see Appendix IB|) 



^{d^^efd,,edi^ = 5^vu{7^) + ^-^^. (37) 



1/2A473 
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Recall once again that the solutions for u and $ repre- 
sent a 'penetrating boundary layer', so one can obtain 
the exact solution with the use of symmetry by taking 
w = at the middle of the domain {y ~ Y/2). For the 
obtaining of the boundary condition \2b\\ and defining TZ 
we evaluate the integral at = 0: 

0(0) 
J {dr,ndfd,fidrj = j sinh^ 9 dO = i(sinh6iocosh6lo-6io). 



Taking into account that E ~ —{T-Voip)y=o and compar- 
ing ([57]) with (pS)) we obtain (the case y ~ Y is similar) 



TZ.^± 







(sinh^o cosh^o - ^o), T^i = ±-p3 ^'o, 

(38) 



TZ = 7^3 + 7^l, 



V. THE FLOW IN A THIN FILM 

In order to describe the flow in a thin film we use the 
simplified version of the equations (fT6l) - ((20| . where we 
accept that the mixture is electroneutral (g = 0) every- 
where but the vicinities of the boundaries. It allows us 
to eliminate from the equations all terms proportional to 
(7, taking them into account only in the boundary condi- 
tions (see Sect. lIVp . The problem describing the averaged 
velocity u = (u, v) and the averaged potential ^ is 

diu + u ■ VoM = — Vop + i^Aou, divo U = 0. (40) 



Ao^ = 0, 



(41) 



where (j4ip corresponds to the continuity equation for an 
electric current in the case of constant conductivity and 
equal diffusion coefficients Dk (see Appendix [C]); ([4T|) 
is not the Poisson-Boltzmann equation that express 
the potential via the charge that was used in Scct. lIVl the 
formal coincidence of these equations should not cause 
misunderstanding. We solve ((40|) , (|4T|) in the rectangular 
domain D ~ {{) ^ x ^ X, ^ y ^ Y} with the boundary 
conditions (^1]). ([Mjl. ^ 



\y=0,Y~ 







(42) 



0,Y' 



TZ^TZiiEo)+TZ3{Eo), (43) 



= 0, 



fa, 



(44) 



Recall that the expression for TZ is given by ((38)) 
and the value of TZ essentially depends on Eq . 

The problem (j41]) . ((44|) has an analytic solution that 
can be presented as Fourier's series. For the further use 
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we give only the following formula (where the sign '+' 
corresponds to y = 0) 



. cxD , , (2k+l)TrY 

The computed graphs of G{x;X,Y) for different values 
of X, Y are shown in Fig.[2l 



Gix;X,Y) 





FIG. 2: The function G{x;X, Y) for X = 1 and the different 
values of Y: (1) F = 1; (2) F = 0.5; (3) Y = 0.2; and (4) 
y = 2 



It is apparent that for the fixed X, Y the sign of ^x\y-Q y 
(and hence the tangential velocity y given by 

depends on the relation between the parameters tpo, Eq. 
For example, for X = Y = 1 and tpo/X = Eq the velocity 
M < on the part of the boundary {0<x<0.1,y = 0}, 
and It > on the rest of it {0.1 < x < l,y = 0}. 



VI. THE NUMERICAL RESULTS 



We solve the Navier-Stokes equations (|40)) with the 
prescribed tangential velocity and the no-leak condition 
at y = 0, Y and the no-slip condition at the rest of 
the boundary (UH), by the employment of the stan- 
dard projection algorithm [H^, [5l| and the finite element 
method. The numerical setting is based on the pack- 



age FreePem-f -I- [52| 
formula (|^5)l for Tp, 



with the use of adaptive grids. The 

„ is not efficient due to its sin- 

v=o, Y 

gularities at x = 0, X] therefore taking into account the 
singularities of derivatives near the vertices we also find 
Ip dH]), (gil) numerically. 

The formulated problem is rather simple, however the 
qualitative properties of its solution strongly depend on 
the relation between the parameters fo, Eq, X , Y. As 
we have already mentioned the direction of the tangen- 
tial velocity on the boundaries y = 0, y is defined by 
(US) (Fig. HI): the velocity is positive on one part of the 
boundary and negative on its remaining part (Fig.[3|) in 



such a way that the particular velocity distribution de- 
pends mainly on the ratio tpo/EQ. It is apparent that 
this tangential velocity causes the rotational motion of a 
large scale. Additional smaller vortices can appear in the 
regions adjacent to the parts of the boundary, where the 
tangential velocity has the opposite sign (Fig. [3]). 




FIG. 3: The sketch of a rotating flow in the flim 

It is instructive to express dimensionless parameters in 
terms of dimensional ones with the use of (|12p , ^3]) , (|38p , 

mi 



T 



- '*^\ F 

Tpl'^l — — Azjiio, 

1 pv* 

^ 135p2j.*3 \rTci) " ' 

1/2 / g*£.2 \ 1/2 



(46) 



£*RT 
2c*gF^ 



V 4RTc*r 



Eo ^ 1.5. 



pa 
ECS' 



Eo = 



'out out 

e* £ 



aE 



We perform our computations for the experimental val- 
ues of parameters for a liquid film motor taken from 
all used values are listed in Tables [11 [Til IIIII It is apparent 
that the velocity {a/T)TZi (that is similar to the classic 
electroosmosis) is significantly less than the tangential 
velocity on the boundary {a/T)TZz that appears due to 
the averaging over the film thickness. Therefore in the 
computations we have not taken TZi into account. 

TABLE I: Dimensional parameters 



Symbol 


Description 


Value 


fo 


difference of potentials 


20 V 


a 


length 


10"^ m 


Eowt 


electric intensity 


30000 V/ m 
lO-^mVs 


V* 


kinematic viscosity 


£,) 


absolute permittivity 


8.85- 10"'^C/(V-m) 


£* 


water permittivity 


78.3 eS 


^out 


air permittivity 


l.OeS 


p 


water density 


10^ kg/m^ 
10-*mol/m^ 


C = c% 


ion concentration 


F 


Faraday constant 


9.65 • 10* C/mol 


R 


universal gas constant 


8.3J/(mol-K) 


T 


absolute temperature 


293 K 
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TABLE II: Characteristic scales 



Symbol 


Description 


Value 


c 

c 
T 
a IT 

\d 
h = Sa 


electric strengtlis scale 
time scale 
velocity scale 
tangent velocity scale 
tangent velocity scale 
Debye's length 
halfheight 


zuuu v/m 
7.8- 10-2 s 
0.128 m/s 
3- 10-2 m/s 
0.5 • 10-'' m/s 
0.95 • 10-" m 
0.29- 10-2 m 



TABLE III: Dimcnsionless parameters 
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0.19 


-1.0 


0.29 


7.8 ■ 10-" 


0.235 


33.42 


1.0 


1.0 


m 


0.19 


-1.0 


0.29 


7.8 ■ 10-" 


0.235 


33.42 


1.0 


0.5 


m 


0.19 


-0.1 


0.29 


7.8 ■ 10-" 


0.235 


33.42 


1.0 


1.0 



One can see that « 0.09; it gives us a sufficient 
ground to treat 8^ as a small parameter and to use (|16p - 

The following figures show the results of computations 
in a square and in a rectangular domain. Fig. [4] shows 
the isolines for the potential ^(x, y) with the step 0.05. 




FIG. 4: The isolines of the potential if{x, y) (left) and ||-(/)(-, || 

Fig.O demonstrates the streamlines of iplx, y, t) with 
the step 0.002 at the instants t = 10 (« 0.78 s) t = 30 
(« 2.34 s). 
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FIG. 5: The streamlines of ipix, y, t) for t = IQ 0.78 s) and 
t = 30 (si 2.34 s) 

The isolines of the velocity field u{x,y,t) at t = 30 
are given in Fig. [51 After t = 30 the flow is practically 



steady; for the additional control of the relaxation to a 
steady state we calculate the mean-square norm \\'ip{-, t)\\ 
(Fig.©. 
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FIG. 6: The isolines of u{x, y, t) (left) and ^(a;, y, t) a.t t = 30 

More detailed discussion of the computational results is 
given in Sect. lVIll Here we just mention that Fig.[5]shows 
the initial appearance of two co-rotating vortices. Later 
on, these two vortices merge into a single vortex that 
represents an almost steady rotating flow in the whole 
domain. For the considered parameters the transition 
(relaxation) to the final steady flow takes around 2 s. 

In addition to the computations in a square domain, we 
perform the computations in rectangular domains with 
different Y . In all cases X > Y we observe a flow struc- 
ture similar to the shown in Fig.O the initial appear- 
ance of two vortices with the subsequent forming of an 
unified steady rotating fiow. For example, the flow for 
a: = 1, y = 0.5 at the instants t = 7 (« 0.546 s) < = 30 
(« 2.34 s) is shown in Fig.[7l 




FIG. 7: The streamlines of i){x,y,t) for X = 1, y = 0.5 at 
t^7 0.546 s) and t ^ 20 1.56 s) 

Fig. [5] shows the flows for the square domain with the 
deliberately smoothed angles (the curvature radius is 
0.1). One can see that the singularities in the electri- 
cal field near the vertices do not alter the flow structure. 
In these computations we keep the boundary conditions 
at X = 0,Ar the same, while on the rest of 
the boundary we introduce physically similar conditions. 
In these computations ^ = on the part [A, B] of the 
boundary, and Ip = (po on [C, D] . The external electric 
field acts in the y-direction. On the rest of the boundary 
the tangential velocity component is proportional to the 
tangential derivative of the potential (similar to ^^)- 
We also keep the no-leak condition valid on the whole 
boundary. 
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FIG. 8: The streamlines for ■ip{x,y,t) at t = 10 (=i 0.78 s) 
and t = 30 (ss 2.34 s) 

We have already mentioned that the tangential veloc- 
ity at the boundary is determined by the relation be- 
tween the parameters (^sq, Eq, X, Y (see (05])) with one 
possible flow regime shown in Fig.[3l In order to con- 
firm its existence we present in Fig. |9] the results for the 
parameters: ^pn = —0.1; Eq = 0.19; X ^ 1; Y ^ I. 
One can see there the isolines of the potential with the 
step 0.01 and the streamlines at t = 200 with the step 
0.0002. The shown flow regime is almost steady: the 
norm ||?A(-,t)|| = 0.001442 in the interval 160 < i < 200 
changes only in the last digit. 
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FIG. 9: The isolines of the potential ifiix, y) (left) and the 
streamlines for 'ip{x, y, t) at t = 200 (si 15.6 s) 

In Fig. [9] the tangential velocity at the boundary y = 
changes its sign at x = Xq « 0.2. The computations 
show that the additional vortices in the angles of the 
domain do not appear if Xq < 0.1. In particular, for 
X ^ 1, Y = 1, Eq = 0.19 the generation of the rotating 
flow takes place when |(^o| > 0.6. It also interesting to 
see the differences between the distributions of potentials 
(c/. Figs-Handd]). 

VII. DISCUSSION 

1. The existence of the discovered EHD rotational flow 
may be expected since it can be generated by the tangen- 
tial velocity at the boundaries. Nevertheless we should 
emphasise once more that a rotational flow appears as 
the result of the applying of constant fields £^out and 
ipo: as it appears in the experiments [H, 01- 



2. The important result of this paper is the ob- 
tained in Scct. lIVh elation between the tangential velocity 
at the boundary and Reynolds stresses. Our averaged 
equations (fTB)) -j^D|) are almost identical to the derived 
in H, d, Ifl HiTalthough we used different boundary 
conditions ©-((Til). In [ii,i,ll3 electrokinetic insta- 
bility for the solutions corresponding to inhomogeneous 
conductivity were studied. The Reynolds stresses terms 
were also derived in these papers, however they had been 
neglected due to their smallness. In our model pB)) - (PDl) 
the situation is right the opposite. Reynolds stresses rep- 
resent the main reason for the appearance of the tangen- 
tial velocity near the boundaries. One can also see in 
Appendix [X] that our averaging method is more detailed 
than the one given in [1, 0, H, H, 0] ■ 

3. A full quantitative comparison of our results with 
the experiments P, Q is impossible, since the key infor- 
mation about the values of some crucial parameters {e.g. 
about the thickness of a film) is absent in these papers. 

4. The qualitative comparison of our results (Figs.[5l 
[3 [5]) with the flow pictures in [l|, Q shows a good agree- 
ment: both in the experiments and in our computations 
one can observe the appearance of the rotational flow, 
growing to its stationary state during the time-interval 
of the order of 2 s. This fact opens the opportunity for a 
fast switching between the directions of a rotation as has 
been proposed in [H, 0] . The magnitudes of rotational ve- 
locities in our results and in the experiments are similar 
(around 3cm/s, at least near the boundaries). 

5. In the experiments the flows with one vortex and 
with two vortices can be observed. Our computations 
show that only one steady vortex can exist. Our results 
show (similar to the experiments) that there are two co- 
rotating vortices in the rectangular film with the ratio 
of sizes 1:2. However our computations also show that 
such a flow is not steady, it finally transforms to the flow 
with single vortex (Fig. [7]). However the authors P,0| do 
not mention whether or not the observed flow with two 
vortices is steady. This contradiction can be resolved 
provided that the experimental observations correspond 
to an unsteady flow. 

6. The experimental rotating flow [l|, 0| appears only 
for some critical values of the electrical field Eq, which 
depend on (pQ . The authors P, 0] mistakenly stated that 
EqLPq = const. Their graph of this function in two loga- 
rithmic scales indeed represents a straight line, however 
its slope is not —1. For our model a rotational 
flow also appears only for the certain values of parame- 
ters. The rough estimation of these parameters follows 
from (see also the comments to Fig.lH]). The rota- 
tional flow with one vortex appears when the tangential 
velocity changes its sign at the point x < O.IX. 

7. The experimental speed of the rotation does not 
depend on the viscosity f* , while the formula for the 
tangential velocity gives {a/T)TZ^ ~ {v*)~^ . However 
for the liquids with different viscosities (the solutions of 
glycerin in water) the thickness of the films also can be 
different, while the velocity is {a/T)TZ3 ~ {v*)-^h'^. We 
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are unable to compare this formula with the experiments, 
since the data on a film thickness in [H, [2] are absent. 

8. In our model p^ - (|20p and in the numerical re- 
sults (Figs.llHS]) the speed of the rotation decreases to- 
wards the center of a film. It looks natural, since the 
cause of this rotation is the tangential velocity at the 
boundary (see the boundary conditions (|^^ ). In con- 
trary, the results [l|, Q show that the speed of the rota- 
tion increases towards the center of a film. On the basis 
of this fact the authors of [J, Q deny electrokinetic ef- 
fects at the film edges as the possible mechanism that 
causes the rotation. However one can propose a num- 
ber of possible explanations for this discrepancy. First, 
it can be the incompleteness of our mathematical model 
that does not consider the surface tension and the de- 
viations of the free surfaces of a film from the planes. 
Second, our mathematical model describes the averaged 
velocity field that differs from the real three-dimensional 
velocity distribution (see (jA17p ). Due to the accepted 
electroneutrality of the mixture (almost everywhere ex- 
cept in the vicinities of the boundaries) the taking into 
account the three-dimensionality of a flow can produce 
the decreasing of the rotation for the layers of a film near 
its boundary. At the same time it is unclear whether 
the data in P, 0] represent the average rotation speed 
or the speed of the rotation of the layer {e.g. the free 
surface) of a film. Third, a more complete mathemati- 
cal model has to consider the Joule heat that naturally 
appears in a weakly conductive liquid under a significant 
electrical current (0.2 mA x 20 V = 4 mW) . The resulting 
nonuniform temperature can cause strong inhomogeneity 
in viscosity and the permittivity of a solution. We should 
recall here that the changing of temperature in the inter- 
val 15-35 °C produces the changing of water permittivity 
Er in the interval 81.9-74.8 {dSr/dT w 0.35). For a strong 
electrical field it can produce a significant pondermotive 
force (l/2)Ve(V(/j)2. 

9. Our model of a rotational flow looks more realistic 
then the heuristic hypothesis of (T], Q on the changing of 
the orientations of water molecular dipoles by an external 
electrical flcld. 

10. The rotating flow in our model is caused by the 
tangential velocity applied at the boundaries. This ve- 
locity has opposite directions at the different parts of 
the boundary. Therefore it is interesting to study more 
systematically the vortex flows that appear at various 
critical values of the applied tangential velocity. 

11. Our model (fTB|) - (P7)) represents only a simplest 
asymptotic model of the flow near the boundary. There 
is a serious potential for the development of this the- 
ory. Here one should keep in mind that the modelling of 
EHD processes in micro-scales represents a rather com- 
plex problem due to the broad spectrum of various phys- 
ical phenomena such as electrokinetic effects (electroos- 
mosis, electrophoresis, etc.), the effects of diffusion, the 
chemical reactions both in a solution and on electrodes, 
the mass-transfer by an electric field, the Joule heat, con- 
vection, Taylor-Aris dispersion, etc. In particular, it is 



unclear weather we can consider the equilibrium Boltz- 
mann concentrations cb ~ 10~^mol/m'^ or we have ions 
of only one sign near the boundaries. 

12. It is especially important to explain the connec- 
tion between our model and the EDL-theories for strong 
external electrical fields [H, HI, [13, [H, Hi, [13, IHI . In 
our model the rotating flow is caused by the edge effects 
at the boundaries j/ = 0, Y , where simplifled bound- 
ary conditions lead to the estimation of the value of 
n = TLiiEo) + TLziEo) (gS]), dM]). At the same time 
this simplified model can be upgraded with the use of 
contemporary EDL-theories (see also the references on 
pp. [2 [4]). This rather complex task can be undertaken if 
the industrial applications of the liquid film motor flows 
appear. Here one can go ahead with the full solution of 
the problem that must include the exact evaluation of 
TZi{Eq) and 'R.j,{Eq) and the correction of the assump- 
tion E ~ const in (see our remark on p. [5]). To 
achieve such a goal one should describe an interface flow 
more precisely, which is possible only with the use of 
EDL-theories. In general, the creation of a full industrial 
level model requires to reconsider or upgrade all results 
of Sect.HVl 

13. In practical applications the liquid film motor flows 
can be used for the micromixing in microfluidic devices. 

14. The general signiflcance of our results for the fur- 
ther developments of microhydrodynamics may consist 
in the revaluation of the role of the considered classical 
effects in the micro- and nano-scale processes. 
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APPENDIX A: THE AVERAGING PROCEDURE 

The averaging of H]), ([S])-®, that takes into account 
the boundary conditions (|5|)- pT|) . gives the exact but not 
closed system of equations 

5^(dtU + u ■ Vqm) + divo('U ® u) = 

= -J^VoP + S^i^Aqu - qVoTp - qVo^, (Al) 

divoM = 0, (A2) 

eAo^ = -q, (A3) 
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dtCk + u ■ VoCk + divo(MCfc) + Dk divo ik = 0, (A4) 

ik = -Dk{VoCk + efc7CfeVo^ + CkjCk'^aip)- (A5) 

In order to obtain the closed system with the precision 
0{6'^) we use the decompositions (fT5|) to calculate the 
terms 

m®m = mO0m" + O((52), (A6) 

W^ = W+'WWWTpW) + o{s^), (A7) 

SfeV;;^ = (^2 + <52Hi)Vo(^" + <52(^i) + 0{S^). (A9) 

For the main terms in (|15p equations ([2]), ([I])-® and 
condition pH]) yield 

(q" + $«R^" = 0, g" + $<' = ^efc(c°+B?.), (AlO) 

fc 

ajo = 0, fk+n^d..4 + eki{cl+4)d^^'- (All) 

(Tl!+7^)L=±i=0. 

Equations (|IlO)) . ((IlT|) give (9^(^ 0, djj! 0, 7° 
cJzC^. It is clear that if dzf = 0, then / = and / = /. 
Hence 

^ = 0, 2^ = 0, $«=0, (A12) 

V'=r. 4=4, q^=t. 

The use of (|A12p transforms the expressions (jATjl - ljAO]) 
to the form 

W^ = OiS'), (A13) 

^ = S^^^ + 0{d^), ^kV^^OiS^). 

From (jA12p we obtain the equation for 

- g°Vo^° = 0, 

which is required for the calculation of (|A6[) with the 
precision 0{5^). In particular it means that we can make 
the replacements q = if + 0{S^), Tp = Tp'^ + 0{5'^) and 
wq, Uq can be foimd from the equations 

1^9,,^° - qVo^ = 0, (A14) 
divo(ll" + M°) + d,w" = 0, (A15) 



with the boundary condition 

w° + w°=0, z = ±l. (A16) 

We assume that vf = 0. The integration of (|A14p - (|A16p 
yields 

u'' = g'{z)U, = -g{z)diYoU, uU^qVo'P, (A17) 

divo«"=0, = i(z3_z), ^=0, 7M=0, 

where we have used the notation (^0)1 for C. 

One can notice that we do not require v!^ to satisfy 
the boundary condition This condition is required 
only for u. The equality 9^1*'^ = at z = ±1 leads to 
U = Q that is not true. In the exact problem one should 
consider a boundary-layer solution at z = ±1 and assume 
the absence of the charge (g = 0) at the boundary. In the 
opposite case the action of a tangential to the boundary 
external field creates the stresses related to Maxwell's 
electromagnetic stress tensor. 

The use of (jAl7[) gives the expression for (|A6P 

fi^ = ^^((7®i7) + 0(<52), y^zj=^. (A18) 

The calculation of uck is based on the next approxima- 
tion for the equations (H])-® 

dt4 + It" • Va4 + w°9,c0 + divo il + dJl = 0, 
il = -Dk(Va4 + ekl4^o^''), 

II = -Dk(dA + e.7(4a.</^° + 4dzV^)) , 

£(Ao/ + 9..^^) = 
or taking in account (jA12p 

dt4 + (m° + • Vo4 + divo 4 + dJl = 0, (A19) 

4 = ^Dk{Vo4 + ekl4^o^), 
ll = -Dk[d,cl + ekl4dz^'), 

The last equation shows that dzz^f^ does not depend on 
z, so to find cj, we obtain the equation 

• Voc° - Dkdzzcl = (A20) 
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with the boundary conditions that follow from (llOp 

For the calculation of uck one can take Ck = + 0{6^) 
in (jA20p . since the required precision for (jASP is 0{S^). 
It allows us to integrate the equation (|A20p 



Finally we obtain 



(A21) 



"fc = {9'{z){go{z) - go{z))^ 



Its substitution into the equation ([16]) for w shows that 
p depends on y only 



It follows that the expression in braces is vanishing and 
(IBTI) gives dSZl). 



APPENDIX C 

Let us consider the case when the values of all diffusion 
coefficients are the same {Dk = D). The multiplying of 
^4p by Ck and combining the results yield 



APPENDIX B 



dq 



A\Ya{uq) — L'divo(Vog + cVqi^) = 0, 



Let us show that in the case ((29|) equations (fT6|) can 
be integrated. The use ([M)) and (PO)) gives the velocity 
component ?I 



Edy (q^dy^) = -(^^a^p + i^S^dyyll - ^-B . (Bl) 



Taking in account that q, <f> and u depend only on y we 
get 

-S^P = ^^Edy {7fdy<^) - v5''dyyU + qE^ X + H {l)) . 



where a = D e^7Cfe is the conductivity of a mixture; 
we have also taken (jA2ip into account. By virtue of (|A2ip 
the electroneutrality ^ = leads to g = everywhere 
except the boundaries. Hence, in the case Ck = cb (see 
([31]) ) we arrive to the equation ([4T|) . 

Recall that the requirement of the equality of all dif- 
fusion coefficients represents a strong restriction. In par- 
ticular, the difference between the diffusion coefficients 
leads to the participation of the term uq (linked to the 
Taylor-Aris dispersion [1, |3j Hi S Q ) i^^to the electroki- 
netic instabilities. 
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